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Preface 


I  chose  this  topic  because  I  am  interested  in  fluid  flow 
but  my  main  field  of  study  at  AFIT  was  structures.  Solving 
a  fluid  flow  problem  using  a  structural  method  satisfied  the 
requirements  of  both  the  school  and  myself.  I  would  like  to 
thank  my  advisor  Capt  J.E.  Marsh  and  iny  typist  Pat  Sawdy  who 
patiently  helped  me  to  finish  this  thesis. 


Dennis  L.  Hunt 
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Abs  L  ra  c  l 


The  finile  element  method  was  used  to  solve  the  two- 
dimensional,  non-linear,  small-disturbance,  velocity-potential 
equation  for  steady  transonic  flow  over  a  thin  airfoil.  Two 
finite  element  upwind  techniques  were  investigated  to  see  if 
either  could  accurately  model  the  supersonic  (hyperbolic) 
zone  that  is  embedded  in  the  subsonic  flow  field.  The  two 
techniques  are:  upwind  functions  and  an  alternative  integra¬ 
tion  scheme.  Both  techniques  used  Galerkin's  Method  of  Weighted 
Residuals,  but  differed  in  the  supersonic  region. 

The  upwind  method  involves  adding  an  upwind  function  to 
the  weight  function  in  order  to  weight  the  upstream  nodes  of 
an  element  more  than  the  downstream  nodes.  The  alternative 
integration  method  involves  Galerkin's  method  for  all  ele¬ 
ments.  In  the  hyperbolic  region,  the  elemental  stiffness 
matrix  is  integrated  only  over  the  area  inside  the  forward 
mach  cones  of  the  elemental  nodes.  Both  these  methods  account 
for  the  physics  involved  in  supersonic  flow.  That  is,  the 
solution  at  a  point  in  supersonic  flow  can  only  be  influenced 
by  points  inside  the  forward  mach  cone  whose  apex  is  located 
at  that  point. 

Neither  of  these  methods  produced  results  that  agree  with 
experimental  data  or  other  solutions.  The  alternative  inte¬ 
gration  method  never  converged.  The  upwind  method  converged, 
but  did  not  converge  to  an  acceptable  solution. 


vi  i 


INVESTIGATION  OF  UEWJN1J  SCHEMES  FOR 


FINITE  ELEMENT  ANALYSIS  OF  TRANSONIC 
FLOW  OVLR  THIN  AIRFOILS 

I.  Introduction 

In  transonic  flow,  nonlinear  equations  with  changing 
characteristics  pose  major  problems  in  solving  the  governing 
equations.  These  equations  need  to  be  solved  since  many  of 
today's  high  performance  aircraft  encounter  some  form  of 
transonic  flow.  Over  the  last  fifteen  years,  investigators 
have  published  many  papers  in  this  area. 

Transonic  flow  calculations  are  important,  because  in 
this  region,  violent  oscillatory  motion  is  often  encountered. 

If  the  effects  of  transonic  flow  are  not  incorporated  into 
an  aircraft's  design,  the  results  could  be  detrimental  to 
its  mission.  Whether  it  be  an  F-16  going  from  subsonic  to 
supersonic  flight,  or  the  tip  of  a  helicopter  blade,  the 
need  for  transonic  flow  calculations  are  evident  in  today's 
high  technology  aircraft. 

Most  of  the  work  done  in  transonic  flow  involves  the 
use  of  potential  flow  theory  and  finite  difference  methods. 

The  velocity  potential  equation  is  used,  because  for  steady, 
irrotational ,  frictionless,  isentropic  flow7  of  a  perfect  gas, 
it  is  a  single  equation  which  satisfies  the  law  of-  conserva¬ 
tion  of  mass,  Newton's  second  principle  of  motion,  and  the 
laws  of  thermodynamics.  Finite  difference  methods  are  used, 
because  they  have  been  around  for  a  long  time,  and  have  proved 
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to  be  a  realistic  wny  to  solve  fluid  flow  problems.  Only 
recently  has  the  finite  element  method  ( II'. ‘-1 )  been  employed 
to  solve  these  types  of  problems.  Initially,  a  method  used 
for  structural  problems,  the  FEM  is  now  bein';  investigated 
to  see  how  well  it  solves  fluid  flow  problems. 

Fi_ni  te  Kit 1  men  t _ Background/! > re v i  chi s  Work 

The  finite  element  method  (FEM)  is  a  numerical  method 
used  to  solve  partial  differential  equations,  and  has  been 
around  for  about  thirty  years.  Prior  to  the  mid  1960’s, 
the  FEM  was  used  to  solve  for  s La  tic  forces  in  structures. 
The  method  was  used,  because  structures  are  nothing  more 
than  an  assemblage  of  finite  parts  connected  at  a  finite 
number  of  points  or  nodes.  Then,  forcing  equilibrium  was 
enough  to  determine  the  loading  at  each  node. 

When  the  mid  1960’s  came,  the  FEM  had  already  proved 
itself  useful  in  solving  structural  problems  using  energy- 
principles  like  the  Rayl iegh-Ri tz  method,  and  was  now  intro¬ 
duced  to  non-s true tura 1  problems  such  as  fluid  flow.  This 
worked  out  very  well  for  simple  problems  in  subsonic  flow, 
and  slowly  has  been  proved  useful  in  more  complicated  non¬ 
linear  problems  such  as  transonic  flow.  Now,  the  FEM  is 
used  to  solve  many  types  of  boundary  and  initial  value 
problems . 

Mathematicians  have  linked  the  FEM  closely  with  varia¬ 
tional  energy  concepts  with  are  found  in  many  fields  of 
structures.  Since  variational  principles  cannot  be  found 
for  every  problem,  the  weighted  residual  method  can  be  used 
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This  scheme  uses  such  methods  ns  col  1  oca  Lion ,  least  squares, 
and  Galerkin's.  In  this  study,  Calerkin's  Method  of  Weighted 
Residuals  will  be  used  for  tlie  analysis. 

Previous  work  in  solving  field  problems  using  the  FHM 
goes  back  to  1965  when  Zienkiewiez  (Kef  1)  published  a  paper 
on  how  to  solve  Laplace  and  Poisson  equations.  A  few  years 
later,  incompressible  flow'  problems  were  solved  using  varia¬ 
tional  methods  and  produced  encouraging  results.  For  tran¬ 
sonic  flow,  the  non-linearity  of  the  equations  caused  conver¬ 
gence  problems  in  an  iterative  solution  scheme.  It  was  not 
until  1975,  when  Shon  and  Habashi  (Ref  2)  were  able  to  get 
converged  solutions  for  mach  numbers  near  critical.  Since 
then,  others  have  obtained  converged  solutions  for  transonic 
flow  using  different  methods.  Wellford  and  Hafez  (Refs  3 
and  4)  used  Galerkin's  method  with  iterative  solution  al¬ 
gorithms  based  on  a  velocity  approximation.  They  were  able 
to  show  convergence  properties  for  freestream  mach  numbers 
well  into  the  transonic  range.  Chan  and  Brashears  (Ref  5) 
used  the  least  squares  method  of  weighted  residuals  to  obtain 
solutions  for  steady  and  unsteady  transonic  flow.  More 
recently,  Marsh  (Ref  6)  used  a  relaxation  technique  applied 
to  the  iterative  non-linear  term  to  get  converged  solutions. 

Object ive 

In  this  study,  the  small-disturbance  potential  equation 
will  be  solved  using  Lw?o  finite  element  techniques  to  see 
if  converged  solutions  can  be  obtained  in  transonic  flow 
over  a  thin  airfoil.  The  two  techniques  to  be  investigated 
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are  an  upwind  method  used  by  Christie  (Ref  7)  and  Heinrich 
(Ref  8),  and  an  alternative  integration  method  suggested  by 
J.  Marsh.  Neither  of  these  methods  have  been  tried  in  con¬ 
junction  with  transonic  potential  flow. 

The  upwind  method  involves  developing  upwind  functions 
that  when  added  to  the  weight  functions  give  the  upstream 
nodes  in  the  supersonic  region  more  influence  than  the  down¬ 
stream  nodes.  This  is  done  to  account  for  the  fact  that  in 
the  supersonic  region,  solutions  at  the  downstream  nodes  of 
an  element  cannot  influence  solutions  at  the  upstream  nodes. 
In  other  words,  the  upwind  function  should  negate  the  in¬ 
fluence  of  the  downstream  nodes  of  an  element  in  the  super¬ 
sonic  region. 

The  alternative  integration  method  involves  more  of  a 
physically  intuitive  approach  to  account  for  the  area  that 
influences  the  solution  at  a  node  in  the  supersonic  part  of 
the  flow.  Knowing  that  only  the  area  inside  a  forward  mach 
cone  can  influence  the  nodal  parameter,  the  finite  element 
equations  were  integrated  only  over  that  cone,  and  not  over 
the  entire  element  as  done  conventionally  in  the  finite 
element  method. 
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II. 


Problem  Description 


In  this  study,  transonic  flow  around  a  thin  symmetric 
non-lifting  airfoil  will  be  determined,  using  potential-flow 
theory.  Solving  the  steady  form  of  the  non-linear  small- 
disturbance  potential  equation  using  finite  element  methods 
creates  a  few  problems.  First,  the  infinite  flow  field  must 
be  reduced  to  a  finite  one  so  it  can  be  discretized  into 
elements.  Second,  the  governing  equation  changes  character 
from  elliptic  to  hyperbolic  as  the  flow  changes  from  sub¬ 
sonic  to  supersonic,  and  third,  the  equation  is  non-linear. 

Flow  Field 

Consider  a  thin  airfoil  in  an  infinite  domain  ft  with 
coordinates  (x,y)  set  up  so  that  the  origin  is  at  the  mid¬ 
chord  (Fig  1).  The  freestream  is  steady  uniform  flow  in 
the  x-direction.  The  domain  ft  extends  from  the  airfoil 
boundary  3fta  to  infinity  as  shown  in  Figure  1. 

In  order  to  simplify  the  problem,  a  few  assumptions 
are  made.  The  infinite  domain  will  be  replaced  by  a  finite 
one  with  the  boundary  Sft^  +  8fta  +  <3^|x|>c/2*  The  airfoil 
is  chosen  to  be  thin,  symmetric,  parabolic  and  non-lifting. 
Thin,  so  only  small  perturbation  velocities  are  present. 
Non-lifting,  so  that  circulation  and  the  Kutta  condition  do 
not  enter  the  problem,  and  symmetric,  so  that  the  problem 
can  be  formulated  in  the  half  space.  The  new  finite  domain 
ft  is  shown  in  Figure  2,  and  the  airfoil  planform  is  given  by 
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Figure  2 


y  ~  -f  (X)  =  ?(i  *4X8) 


(1) 


where  t  is  the  thickness  ratio  t/c  shown  is  Figure  2. 

Governing  Differential  Equation 

Assuming  steady,  irrotat iona 1 ,  inviscid,  isentropic 
flow  of  a  perfect  gas  permits  the  use  of  potential  theory 
to  describe  the  velocity  changes  in  the  freestream.  The 
presence  of  a  thin  airfoil  permits  the  use  of  small-distur¬ 
bance  theory.  The  non-dimensional  small-disturbance  poten¬ 
tial  equation  (Ref  9)  in  terms  of  the  non-dimensional  velocity 
potential  is 

^  =  0  (2) 

where  is  the  freestream  mach  number  and  y  is  the  ratio 
of  specific  heats  (y=1.4  for  air).  Equation  2  is  valid  for 
all  points  (x,y)  in  ft. 

Equation  2  is  used  for  the  case  of  transonic  flow  where 
the  non-linear  term  <{>  <j>  becomes  significant  and  is  re- 

j  X  j  XX 

quired  to  describe  the  flow  phenomena  that  occurs.  For  low 
mach  numbers  (i.e.  Mro<0.5)  this  term  is  small  compared  to 
the  others  and  is  often  neglected.  For  the  incompressible 
case  (>1^=0. 0)  equation  2  reduces  to  the  Laplace  equation. 

In  this  study,  the  form  of  equation  2  will  be  used  for  all 
mach  numbers  from  zero  to  just  less  than  one. 

From  partial  differential  equation  theory  (Ref  10), 
eq  2  is  referred  to  as  a  second-order,  non-linear  partial 
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differential  equation  of  mixed  character.  The  equation  is 
second-order,  because  its  highest  derivative  is  a  second 
derivative.  Non-linear,  because  of  the  C  v'  term,  and 

)  X  j  XX 

mixed,  because  the  equation  changes  character  from  elliptic 
to  hyperbolic  as  the  local  velocity  goes  from  subsonic  to 
supersonic,  respectively.  It  is  this  change,  that  constitutes 
the  name  transonic.  This  change  happens  at  some  critical 
value  of  freestream  mach  number  >1  when  the  flow  at  a  point 
near  the  center  of  the  airfoil  becomes  sonic.  As  M  is 

cu 

increased  further,  a  region  of  supersonic  flow  forms  over  the 
airfoil  and  is  called  a  supersonic  bubble.  When  this  region 
gets  big  enough,  a  weak  compression  shock  forms  in  the  down¬ 
stream  part  of  the  bubble  to  allow  the  flow  to  return  to 
subsonic  speed.  Along  the  upstream  boundary  of  the  super¬ 
sonic  bubble  the  local  mach  number  is  equal  to  one  and  the 
coefficient  of  the  <t>  term  in  eq  2  is  zero;  therefore,  the 

j  XX 

equation  is  parabolic.  The  development  for  transonic  flow 
over  ar:  airfoil  is  shown  in  Figure  3. 


Boundary  Conditions 

The  boundary  conditions  for  the  flow  field  come  from  the 
fact  that  the  perturbation  velocities  (u,v)  must  go  to  zero 
at  infinity,  and  there  cannot  be  any  flow  through  the  boundary 
of  the  airfoil  3fta.  In  equation  form,  the  boundary  conditions 
are 


1)  V<)>  — *•  0  AST  — ►  eo 

2)  FOR  (X,Y)  IN  SDa 


(3) 
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Figure  3.  Effect  of  Freestream  Mach  Number  on  Local 
Mach  Number  for  Thin  Airfoils 
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where  r  is  any  distance  from  the  origin,  i  is  the  full 
potential  and  n  is  a  unit  vector  normal  to  the  surface  of 
the  airfoil.  Further  development  of  these  boundary  condi¬ 
tions  is  given  in  Section  III.  Along  the  far  field  boundary 
an  asymptotic  solution  developed  by  Klunker  (Ref  11)  will 
be  imposed. 


III.  Analysis 


In  Section  11,  the  governing  equation,  the  boundary 
conditions  and  the  difficulties  associated  with  transonic 
flow  calculations  were  discussed.  This  chapter  describes 
the  formulation  techniques  necessary  to  solve  the  problem 
numerically  using  the  finite  element  method.  The  reader  is 
urged  to  refer  to  the  appendices  for  additional  information. 

Flow  Field  Discretization 

The  infinite  domain  ft  discussed  in  Section  II  must  be 
replaced  by  a  finite  flow  field  so  that  it  can  be  discretized 
into  a  finite  number  of  elements.  This  can  be  done  two  ways. 
First,  the  far  field  boundary  ft^  could  be  taken  to  be  very 
large,  and  the  actual  boundary  condition,  given  by  eq  3,  be 
imposed.  This  choice  may  not  be  a  good  one,  because  a  large 
flow  field  requires  a  greater  number  of  elements,  which  means 
higher  computational  costs.  Also,  because  the  boundary  con¬ 
ditions  are  of  the  Neumann  type,  the  solution  can  only  be 
determined  to  within  an  arbitrary  constant  unless  a  specific 
value  of  <{)  is  given  for  some  arbitrary  point  (Ref  6). 

An  alternative  technique,  used  in  this  study,  is  to  use 
the  far  field  solution  developed  by  Klunker  (Ref  10).  For 
this  technique,  $  =  <+> ^,p  is  specified  along  the  far  field 
boundary.  Klunker's  solution  satisfies  the  actual  boundary 
condition  (eq  3)  and  is  valid  only  at  points  in  the  far 
field  of  ft.  This  approach  has  been  successfully  used  by 
researchers  using  either  finite  difference  or  finite  element 
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(Refs  2,5,6)  methods.  Klunker's  method  allows  the  use  of  a 
much  smaller  domain,  as  compared  to  the  first  method  dis¬ 
cussed.  This  small  domain  means  fewer  degrees-of-f reedom 
and  lower  computer  costs.  Klunker's  equation  for  the  far 
field  potential  pp  for  a  non-lifting  airfoil,  ignoring 
higher  order  terms  is 

(b  -  J - -  fflpd?  (4) 

Tff  rf/Q  X2H/9Y)  J  7 

c 

where 

and  f(x)  is  given  by  eq  1.  Evaluating  the  integral  gives 

A  =2  .  JCC  .  _ ! _  (5) 

%  3  rr,e  xJ  + y!(i-mS,> 

which  is  a  function  of  the  location  of  the  far  field  points. 

Working  with  a  symmetric  airfoil  at  zero  angle  of  attack 
allows  the  problem  to  be  set  up  in  the  half  space.  The  flow 
field  can  be  designated  by  two  parameters,  Xmax  and  Ymax,  as 
shown  in  Figure  4.  In  order  to  find  optimum  values  for  Xmax 
and  Ymax,  results  of  tests  done  by  Marsh  (Ref  6)  will  be 
used.  Marsh,  solving  this  same  problem  by  a  different  tech¬ 
nique,  varied  the  number  of  elements  in  the  mesh  keeping  the 
element  size  the  same.  He  found  that  the  pressure  distribu¬ 
tion  converged  to  the  assumed  solution,  for  mach  numbers  from 
0.0  to  0.8,  when  the  value  for  Xmax  and  Ymax  were  1.5  chord 
lengths  or  greater.  In  this  study,  the  minimum  dimension 
(1.5c)  will  be  used  for  both  Xmax  and  Ymax. 
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Figure  4.  Flow  Field  Parameters 
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Working  in  the  half  space,  the  domain  of  the  problem 
can  be  broken  down  into  three  regions  as  shown  in  Figure  4. 
Regions  1  and  ill  are  upwind  and  downwind  of  the  airfoil 
respectively  and  can  be  discretized  into  fairly  large  ele¬ 
ments.  Region  11,  over  the  airfoil  will  be  discreti/.ed 
finer,  because  it  is  here  that  perturbation  velocities  change 
the  most.  Over  the  airfoil  is  also  where  Lho  pressure  dis¬ 
tribution  is  required,  so  more  nodes  are  needed  to  accurately 
depict  the  pressure  distribution.  Using  a  thin  airfoil,  the 
boundary  terms  of  the  finite  element  equations  (given  later 
in  this  section)  are  evaluated  along  the  chord  (y=0)  of  the 
airfoil.  Therefore,  the  elements  in  Region  II  extend  down 
to  the  x-axis  and  do  not  terminate  at  the  airfoil  contour. 

The  mesh  shown  in  Figure  5  is  representative  of  the 
meshes  used  in  this  study.  Variations  include  changing 
the  number  of  elements  over  the  airfoil,  and  in  Regions  I 
and  III. 

Finite  Element  Solution 

When  solving  partial  differential  equations  using 
finite  element  methods,  it  is  common  practice  to  transform 
the  governing  equation  into  the  form  of  the  matrix  equation 

[K]HHF]  (6, 

where  (K)  is  the  stiffness  matrix,  {$}  is  the  solution  vector 
and  {F}  is  the  forcing  terra.  The  method  of  weighted  residuals 
will  be  used  to  accomplish  this. 
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The  steady  form  of  the  non-dimensional  sma  1  1  -d  i  s  l  rtibance 
potential  equation  (Kef  9)  for  transonic  flow,  g  i  veil  by  eq  2, 
in  two  dimensions  is 

+  (2) 

for  all  points  (x,v)  in  the  domain  U.  Relationships  between 
the  non-dimensional  parameters  (d,x,y)  and  the  physical  ones 
( 4> , x , y )  are:  x  =  x/c,  y  =  y/c,  and  4>  =  $/uwc,  where  is 
the  freestrcam  velocity  and  c  is  the  chord  length.  If  eq  2 
is  formulated  by  finite  element  methods,  a  set  of  second- 
order,  non-linear  algebraic  equations  will  result.  This  can 
be  simplified  two  ways,  one  using  an  iterative  scheme  to 
linearise  the  non-linear  term;  and  two,  integrating  by  parts 
to  lower  the  number  of  continuous  derivatives  repaired  in  the 
assumed  solution  so  linear  elements  can  be  used.  Both  of 
these  simplifications  will  be  employed  in  this  study,  and  arc 
elaborated  on  in  the  next  few  paragraphs. 

Rewriting  eq  2  as 

[(l-M»)<j>x-Mi(l+yl/2  <£],„+  <j>83=0  (7) 

2 

leaves  0  as  the  non-linear  term.  Using  the  iterative 
> x 

approximation  (Ref  6) 

£  =&  C  .  (8) 

where  superscript  n  denotes  the  iteration,  allows  the  poten¬ 
tial  <t>  to  be  replaced  by  a  sequence  of  potentials  { $ ^ ,  d ^  ,  .  .  .  , 
£n,fn+^},  which  converge  when 


17 


(9) 


where  e°  is  some  small  error  chosen  as  Lhe  convergence  re¬ 
quirement.  In  this  study,  e°  is  chosen  to  be  five  hundredths 
of  a  percent  (Ref  6).  Now,  eq  7  in  iterative  form  becomes 


[(i-mSJC'-m-—-’  C  CL+  C  =  0  (10) 


where  n+1  denotes  the  variable  $  to  be  solved  for,  and  n 
denotes  the  variable  4>  calculated  from  the  previous  iteration 
The  method  of  weighted  residuals  is  a  technique  used 
to  generate  finite  element  equations  when  variational  func¬ 
tions  are  not  available.  This  method  assumes  an  approximate 
solution  of  the  form 


<j>(X>Y)  g  Nl(X,Y)  cj>b  (ll) 

where  the  N^'s  are  functions  that  satisfy  boundary  conditions 
and  the  c^'s  are  the  solutions  at  global  node  points  i. 
Substituting  this  approximate  solution  into  eq  10  results 
in  the  equation  not  being  equal  to  zero,  but  now  being 
equal  to  some  error  e.  As  the  error  approaches  zero,  the 
approximate  solution  approaches  the  exact  one.  Therefore, 
the  object  is  to  get  this  error  as  small  as  possible.  This 
is  done  in  a  weighted  average  sense  by  substituting  e  into 

f &W;.da  =  0  (12) 
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where  is  a  weight  or  test  function.  Substituting  for  c 
in  eq  12  yields 


l  +  ff 
2 


-Win  ^ 

),x+^  )widn=o 


(13) 


which  is  still  second-order  and  requires  to  be  continuous 
through  its  first  derivative  (C^).  Integrating  eq  13  by 
parts  eliminates  the  continuity  by  shifting  a  derivative 
from  $  to  W ^  and  thus  requires  N.j(x,y)  and  W-(x,y)  to  only 
be  continuous  functions  (C^1).  The  resulting  equation  after 
integration  by  parts  is 


-  Ml  1-Mi-Mi 


!«■ 

“XI 


<P,J‘R*  wi,x+  L  wS)d0 


+  i;v4dS=0 


(14) 


where  3ft  denotes  the  boundary  and  n  =  (nx,n^)  is  the  unit 
normal  vector  to  the  boundary  surface. 

When  evaluating  the  boundary  term  of  eq  14,  the  integral 
must  be  broken  down  into  a  series  of  integrals  over  each  of 
the  six  boundaries  shown  in  Figure  6.  Before  evaluating 
these  integrals,  the  boundary  conditions,  given  by  eq  3,  for 
the  governing  differential  equation  must  be  put  into  a  work¬ 
able  form.  Recalling,  the  boundary  conditions 


0  as  r  — ►  CO 


(15) 


and 

v£  •  n*=  o 


along  airfoil  (16) 
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where  the  full  potential  £  is  defined  as 

£=Uoo(X  +  <j>)  (17) 

and  is  the  unit  normal  vector  on  the  airfoil  given  by 

-  _  VF  Cf.x  -!> 

‘  ~  II V F]|  -  IIVFII  <18) 

where 

F  =  f(X)~  Y  —  0 

is  the  equation  of  the  airfoil  given  by  eq  1.  Substituting 
eqs  17  and  18  into  eq  16  results  in  the  boundary  condition 

ir  (|  +  V  <19> 

along  the  airfoil.  From  symmetry  and  steady  flow  in  the 
x-direction  another  boundary  condition 

4>  ~  0  (20) 

•  >s 

exists  for  y=0  where  |x|  >  1/2. 

In  evaluating  the  boundary  term  of  eq  14  for  far  field 
boundaries  1,  2,  and  3,  in  Figure  6,  the  infinity  boundary 
condition,  given  by  eq  15,  should  be  used.  Since  Klunker's 
boundary  condition  q  =  given  by  eq  5  is  used  instead, 

then  VT  will  be  taken  to  be  zero  there.  The  boundary  terms 
due  to  segments  4  and  6,  in  Figure  6,  are  also  zero  when 
evaluated  using  n^  =  0  and  eq  20.  The  only  term  that  docs 
not  go  to  zero  in  the  boundary  term  is  for  segment  5.  Sub¬ 
stituting  eqs  18  and  19  into  the  boundary  term  of  eq  14 
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results  in 


-  Mw^-W^l dX 

*aa  w=o 


Since  $  and  will  be  continuous  functions  along  inter¬ 
element  boundaries,  it  is  possible  to  integrate  eq  14  piece- 
wise  over  the  domain.  Within  each  element,  eq  11  for  a  four- 
noded  rectangle,  becomes 

“tn+t  4  .  rv-rt 

<j>  (X,Y)  =  ]>  Nj ( X  Y )  ® 

1  3=1  J 

or  (22) 

*Tn  4  ,  n 

9  (x,y)-^nk(xy)9k 

where  the  Nj 1 s  and  N^'s  are  now  elemental  shape  functions. 

The  shape  function  for  a  four-noded  rectangular  element  is 
given  in  Appendix  A.  Substituting  eqs  22  and  21  into  eq  14, 
produce  the  finite  element  equations  in  elemental  form.  The 
equations,  in  the  form  described  by  eq  6,  are 

[( I-  Ml)  A,j  +  B;j  -Ml  ~  C4tf )  +  Ml  Dy 

(23) 

+  Ml^~El.J(f>]  cjT'  =  Fl 


where 


AM  =  jjW,XNi.,dXdY 

sx 

B«=II  W^N^dXdY 


22 


Cy(f  >  =  <t> 


dXdY 


'a. 


Dtj  = 


=  (  N^WifJdX 
s=0 


(24) 


E4^1=  <t>"f  N^Ni^Wj 

ydn« 


>x 
y=o 


dX 


»X 

an*  yso 


dX 


where  3fta  denotes  that  these  terms  apply  only  to  elements  on 
the  boundary  of  the  airfoil.  These  matrices  are  evaluated 
for  the  four-noded  rectangular  element  in  Appendices  A,  B, 
and  C . 

All  of  the  parameters  in  eq  23  are  specified  except 
the  weight  function  W^.  The  reason  for  this,  is  that  the 
weight  function  in  this  study,  will  change  depending  on 
the  method  used  and  whether  an  element  is  contained  within 
the  supersonic  region.  For  all  the  upwind  method  elements 
in  the  subsonic  region  (elliptical),  and  for  all  the  elements 
in  the  alternative  integration  method,  Calerkin's  method  will 
be  used.  This  means  that  for  those  cases,  the  weight  function 
will  be  the  same  as  the  shape  function  N^.  In  the  upwind 
method,  the  weight  function  for  supersonic  (hyperbolic) 
elements  will  be 


23 


wv  -  n^+ocu. 


(25) 


where  IK  is  an  upwind  function  and  a  is  a  test  coefficient. 
Both  the  upwind  function  method  and  the  alternative  integra¬ 
tion  method  will  be  discussed  in  detail  later  in  this  section. 

Transonic  Problems 

Chapter  II  explained  what  happens  in  the  flow  as  tran¬ 
sonic  speeds  are  reached.  This  section  explains  how  the 
phenomena  will  affect  the  finite  element  solutions. 

Researchers  studying  the  finite  element  (Refs  12,  13) 
and  finite  difference  methods  for  transonic  flow  have  re¬ 
ported  convergence  difficulties.  Some  believed  the  problem 
was  in  the  small-disturbance  potential  equation  (Rut  6),  but 
Akay  (Ref  12)  reported  convergence  difficulties  for  the  total 
potential  equation.  Others  believed  the  problem  w as  in  the 
Galerkin  formulation,  because  it  does  not  account  for  the  area 
of  influence  of  supersonic  nodes  (Ref  5). 

In  the  supersonic  region,  the  governing  equation  is 
hyperbolic.  From  the  aerodynamics  of  supersonic  flow,  a 
point  in  the  flow  cannot  propagate  waves  forward;  it  can 
only  be  influenced  by  points  that  lie  in  a  region  inside 
the  mach  cone  propagating  forward,  and  can  only  influence 
points  that  are  contained  in  the  downstream  cone.  *  The 
mach  cone  is  defined  by  the  characteristic  curves  of  the 
hyperbolic  equation.  It  is  believed  that  neglecting  these 
phenomena  leads  to  convergence  difficulties  using  Galerkin' s 
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methoJ.  Far  Galerkin's  method,  the  nodes  of  an  element  are 
weighted  the  same  when  the  upstream  nodes  should  have  more 
influence  in  supersonic  elements. 

In  order  to  account  for  the  hyperbolic  behavior,  re¬ 
searchers  have  tried  various  methods.  In  finite  difference 
methods,  special  difference  operators  (backward  differences) 
have  been  developed  to  insure  convergence  of  the  solution. 

In  finite  elements,  there  has  not  been  a  special  formulation 
developed.  Chan  (Ref  13)  changed  the  stiffness  matrix  for 
supersonic  elements  during  the  assembly  process  by  zeroing 
out  the  rows  in  the  stiffness  matrix  corresponding  to  nodes 
giving  downwind  influence.  Chan's  method  gave  converged 
results,  but  only  for  the  least  squares  formulation.  Marsh 

(Ref  5)  modified  the  non-linear  term  4>  n+^  by  putting 

,  x  ,  x 

in  a  relaxation  term.  His  new  term  replaced  ^  vn4>  n+^  by 
R <f>  n<t)  n+^  +  ( 1 — R ) n<^  nU,  where  R  was  the  relaxation 

>  x  >  X  ,  X  ,  X 

coefficient  which  ranged  from  0  to  1 ,  and  U  was  an  upwinding 
factor  which  also  ranged  from  0  to  1 .  The  second  term  was 
then  added  to  the  forcing  vector  in  the  elemental  equations. 
Marsh  got  converged  solutions  for  mach  numbers  deep  in  the 
transonic  range.  His  method  worked  quite  well,  the  only 
difficulty  encountered  was  selecting  the  values  of  R  and  U 
for  any  given  freestream  mach  number.  In  this  study,  two 
new  methods  will  be  tried.  One  used  by  Christie  (Ref  7),  an 
upwind  method,  and  the  other  suggested  by  Marsh,  an  alterna¬ 
tive  integration  method. 
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Upwind  Function  Method 

The  purpose  of  the  upwind  function  is  to  weight  the 
upstream  nodes  of  supersonic  elements  more  than  the  down¬ 
stream  nodes.  The  reason  for  this  was  explained  previously 
in  this  section.  The  upwind  function  idea  was  taken  from 
papers  by  Christie  (Ref  7)  and  Heinrich  (Ref  8).  Christie 
applied  upwind  functions  in  second-order  equations  with 
significant  first  derivative  terms,  and  Heinrich  used  it 
for  the  convective  transport  equation.  Neither  one  of  these 
parallels  the  transonic  potential  flow  problem  in  this  study, 
but  the  method  seemed  worthwhile  to  investigate. 

Using  Galerkin's  method,  the  weight  functions  W-  are 
assumed  to  equal  the  shape  functions  N • .  Convergence 
problems  arise  when  transonic  flow  is  present,  because  the 
shape  functions  weight  each  node  in  an  element  the  same. 

For  elliptic  elements  (subsonic)  Galerkin’s  methos  (W^-N^) 
is  used,  but  for  hyperbolic  and  transition  elements  (super¬ 
sonic)  an  upwind  function  is  added  to  give 


Wi=Ni+  auv 


(hyperbolic) 


Wy  =  Ny.  +  £  U  v 


( transition) 


where  U^  is  the  upwind  function  and  a  is  a  test  coefficient. 
Transition  elements  are  ones  that  contain  both  elliptic  and 
hyperbolic  nodes. 


Shape  Functions.  Elemental  shape  functions  (Ref  1), 
must  be  equal  to  one  at  node  i  and  aero  at  the  other  three 
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nodes.  The  shape  functions  for  the  bilinear  rectangle  used 
in  this  study  are 

Ni=5'(l+tlt)(l+?,?)  (26) 

where  (£,n)  are  the  variables  in  local  coordinates,  and 
( ^  , n  ^ )  are  the  local  coordinates  of  node  i,  i  =  1,2, 3, 4. 
Figure  7  shows  the  relationship  of  shape  functions  within 
the  element.  By  superposition,  the  nodal  shape  functions, 
eq  26,  are  added  together  to  give  the  curve  of  the  elemental 
shape  function  in  Figure  3.  Notice  that  the  curve  is  con¬ 
stant  over  the  element  boundary.  This  shows  that  all  the 
nodes  are  weighted  the  same. 


Upwind  Function.  The  purpose  of  the  upwind  function  is 
to  weight  the  upstream  part  of  an  element  in  the  supersonic 
region  more  than  the  downstream  part.  This  will  be  done  by 
making  a  piecewise  parabolic  function  of  the  form  (see 
Fig  9) 


u(= 


[Cf-H,)'+K,](l+*?i'7)  -1<?<0 


(27) 


^[(t-H,)2  +  K2](i+m'!)  o<r<i 


where  (H,K)  are  the  coordinates  of  the  vertex  of  the  parabola 
and  P  is  half  of  the  latus  rectum  (twice  the  distance  from 
the  vertex  to  the  focus).  Figure  10  shows  a  plot  of  = 

Nf  +  alT^  (u=l)  so  the  effect  of  the  upwind  function  can  be 
seen.  Notice,  now  the  upstream  part  of  the  element  is 
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weighted  more  than  the  downstream  part. 


Elemental  Unwind  Equations.  Substituting  W-  =  N-  +  uU • 
- - - i -  °  1  1  i 

into  the  elemental  equation,  given  by  eq  23,  yields  a  net; 
set  oi  equations  to  be  used  for  hyperbolic  elements: 


[(l-^)(A:.rHXAUlj)+Bq+aBUlj-M^L|I(Cu^aCUl) 


+  ( Dq+ecDUij)  +  M Eq  +CCE UkJ) ] ]  =  - ( c  -V  cc FUL ) 


(28) 


where  A-  •  ,  B--,  C--,  D  •  .•  ,  E-  .  ,  and  f.  were  defined  by  eqs  24 

^  J  ^  J  ^  J  ■*- 

and 


AUii~{{N  J  xUi>x  dxdvj 


n 


bum~  J[Ni„U^d*dvi 

n 

CUii(f)=CjfN.„Nj,K  Ui.,dxdy 

.Cl 


DU 


irR.ua, 


dx 


r-o 


EU,i<f)=  €  \  N^Nlf,  Utf, 

dfl* 

FUr  J  ui  L 


S  ~-o 


dx 


bXt  *.  ^-0 


(29) 
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where  again  means  this  integral  is  only  applied  for 

elements  on  the  airfoil.  Appendix  B  evaluates  eqs  29  in 
matrix  form  for  a  bilinear  rectangular  element. 

For  the  upwind  method,  a=0  will  be  used  for  elliptic 
elements,  and  a  greater  than  zero  for  hyperbolic  elements. 

The  elements  that  have  some  nodes  elliptic  and  some  hyper¬ 
bolic  (transition  elements)  will  use  a  value  of  3  greater 
than  or  equal  to  zero.  Note  that  a  =  3  =  0  reverts  the 
elemental  upwind  equations  (eq  28)  t:o  the  elemental  equations 
given  by  eq  23. 

Alternative  Integration  Method 

This  method  is  based  upon  physical  intuition  and  uses 
Galerkin's  method  for  all  elements.  A  modification  is  made 
for  hyperbolic  elements;  the  integrals  (eq  24)  are  integrated 
only  over  the  element  area  contained  inside  the.  forward  mach 
cone  (i.e.  not  over  the  entire  element  as  done  in  Appendix  A). 

To  find  the  form  of  the  elemental  stiffness  matrix, 
the  equations  of  the  mach  lines  must  be  found.  The  mach 
lines  are  given  by  the  characteristic  lines  of  the  governing 
differential  equation  given  by  eq  2  as 

^  +  ^4>„=  o  (30) 

where  1/c^  =  1  -  M^( 1  +  y)b  n.  Since  u  =  $  can  be 

calculated  from  the  previous  iteration  and  the  element  is 

2 

small  compared  to  the  flow  field,  1/c  will  be  assumed 
constant.  When  the  e  ,  term  is  negative,  eq  30  becomes 
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(31) 


<J)  —  rk  c|)  —  0 
‘iVi  C2 

2 

which  is  hyperbolic.  With  1/c  constant,  the  characteristic 
lines  (Ref  10)  are  defined  as 

Y  =  CX-f-  d  ;  Y--CX  +  d  (32) 

which  must  now  be  transformed  to  local  coordinates  by  the 
transformation  equations 


x-xc 

a 


and 


Y~YC 

b 


where  xc,  yc,  a,  and  b  are  shown  in  Figure  7.  The  charac¬ 
teristic  (mach)  lines  now  take  the  form 


where 


9=M'f  +  B  ;  +  B 


M=C-|- 


B  -  C 


Xr. 

b 


b 


(33) 


Now  that  the  equations  for  the  mach  cone  are  known  for  any 
point  in  the  supersonic  region,  the  influence  on  the  stiff¬ 
ness  matrix  ^k"]  can  be  determined.  is  defined  (Ref  1) 

as  the  force  at  node  i  due  to  the  unit  potential  at  node  j. 
Therefore,  only  terms  in  which  node  i  falls  into  the  mach 
cone  at  point  j  have  influence  in  the  solution. 
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When  an  element  aspect  ratio  b/a  is  greater  than  or 
equal  to  the  slope  c  of  the  mach  lines,  the  elements  take 
the  form  shown  in  Figure  11a  and  its  stillness  matrix  Lakes 
the  form 


Similarly,  when  b/a 
stiffness  matrix  is 


11 

0 

0 

0 

21 

0 

0 

K24 

31 

0 

0 

K34 

0 

0 

K44 

less 

than 

1  C, 

see 

11 

21 

0 

0 

0 

0 

0 

0 

0 

0 

K34 

0 

0 

K44 

where  K—  in  both  cases  is  calculated  for  a  bilinear  rec¬ 
tangular  element  in  Appendix  C. 

Transition  elements  pose  more  of  a  problem,  because  a 
linear  interpolation  scheme  must  be  imposed  to  find  the 
value  of  n*  where  the  element  changes  from  hyperbolic  to 
elliptic.  If  the  slope  of  the  characteristic  curve  passes 
through  £  =  -1  below  n*,  that  intercept  must  be  used  as  the 
limit  of  integration  in  the  n-direction,  otherwise,  n*  will 
be  used.  In  either  case  the  element  shown  in  Figure  11c 
has  a  stiffness  matrix  of  the  form 
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a)  Supersonic  Element,  c  >  a/b 


b)  Supersonic  Element  c  <  a/b 


c)  Transition  Element 


F igure  1 1 . 


Mach  Cones  at  Nodes  of  Supersonic  and 
Transition  Elements 
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Again,  Appendix  C  calculates  the  values  for  K^j . 

Summarizing  the  alternative  integration  method,  the 
elliptic  elements  use  Galerkin's  method  integrated  over  the 
entire  element.  The  hyperbolic  and  transition  elements  are 
integrated  over  the  area  inside  the  forward  rnach  cone.  For 
example,  for  b/a<c,  the  limits  of  integration  for  K^,  for 
the  area  shown  in  Figure  12,  are 

I  1 

j  |  F(T,9>dr?dt 

-I  M,1*B 
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IV  Results 


Transonic  flow  over  a  six  percent  thick  airfoil  was 
solved  for  using  the  two  upwind  techniques  described  in 
Chapter  III.  For  each  technique,  different  parameters 
were  varied  in  order  to  study  the  convergence  behavior 
of  the  solution.  Before  the  transonic  cases  were  tried, 
the  incompressible  (>1^=0)  case  and  subsonic  cases  were 
tried  to  see  if  the  finite  element  equations  defined  in 
Appendix  A  accurately  modeled  the  flow. 

For  incompressible  flow,  the  solution  shown  in  Figure 
13  was  found  by  directly  solving  the  finite  element  equa¬ 
tion.  This  solution  is  compared  with  the  exact  solution 
(Fig  13)  for  two  discretizations.  Grid  A  was  the  same  as 
shown  in  Figure  5  and  Grid  B  was  similar,  except  it  had 
ten  divisions  over  the  airfoil  instead  of  eight.  Taking 
the  value  of  the  coefficient  of  pressure  in  an  element  to 
be  at  the  mid  point  between  nodes,  Figure  14  shows  that 
the  finite  element  solution  "sely  approximates  the  exact 
solution . 

Subsonic  cases  were  tried  for  freestream  mach  numbers 
between  0.2  and  0.8.  All  these  cases  were  below  the  criti¬ 
cal  mach  number  M  ,  and  converged  within  four  iterations. 
They  also  agree  with  linear  theory  until  M  got  hear  0.8, 
when  the  non-linear  term  started  to  show  its  effect.  Figure 
15  shows  the  pressure  distribution  for  >1=0. 2,  0.4,  0.6,  0.8, 
and  Figure  16  shows  the  difference  between  the  linear  (Ref  6) 


Distance  from  Leading  Edge 

-  Exact  Linear  (Ref  6) 

O  FEM  (Grid  B) 


Figure  16. 


Comparison  of  Finite  Element  and  Exact 
Linear  Solutions  for  M  =0.8 


and  finite  element  solutions  for  Mw=0.8.  Notice  that  the 
difference  between  solutions  is  slight,  even  at  11^=0. 8. 

Using  Galerkin's  method,  without  upwind  techniques, 
mach  numbers  above  0.8  were  run  to  find  the  critical  value 
of  Mw.  Figure  17  shows  the  pressure  distributions  found 
for  the  mach  numbers  that  converged.  The  critical  free- 
stream  mach  number  was  found  to  be  between  Moo=0. 83  and  0.84 
for  a  six  percent  thick  parabolic  airfoil.  Notice  that  the 
solution  does  converge  for  MOT=0.84  and  0.85.  This  is  be¬ 
cause  the  supersonic  region  is  small.  Figure  18a, b  shows 
the  transition  and  hyperbolic  element  in  the  flow  field 
for  M^O.84  and  0.85. 

In  Figures  19  and  20,  convergence  behavior  is  shown 
using  Galerkin's  method  for  >1^=0. 85  and  >1^=0. 86.  At  11^=0.85 
the  solution  converged  after  seven  iterations,  and  at 
11^=0.86,  the  solution  diverged  because  of  the  presence  of 
two  hyperbolic  and  four  transition  elements.  As  the  free- 
stream  mach  number  increased  from  zero  to  0.85,  the  number 
of  iterations  required  to  get  a  converged  solution  increased, 
and  above  M^O.85,  the  iterative  scheme  diverged.  This  be¬ 
havior  is  graphed  in  Figure  21  for  convergence  criteria  in 
the  sense  of  eq  9. 

These  calculations  proved  that  the  finite  element  method 
as  formulated  worked  until  transonic  flow  developed.  The 
techniques  discussed  in  Chapter  III  were  then  tried  to  see 
if  convergent  solutions  occur  for  M  <  Mm  <  1.0  which  is  in 
the  transonic  region.  The  rest  of  this  section  discusses  the 
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Convergent  Results  for  M,o=.82, 
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Figure  20.  I  Lora  Live  Divergence  Behavior  for  -  0.S6 

47 


results  obtained  for  the  upwind  and  alternative  integration 
methods . 


Upwind  Method 

The  upwind  function  method  discussed  in  Chapter  III  was 
tried  for  several  variations  of  the  test  coefficients  a  and 
8,  the  area  under  the  upwind  function  curve,  the  freestream 
mach  number  Mm,  and  the  y-direction  upwind  influence  (BU— ). 


Y-Direction  Influence.  The  upwind  functions  given  by 
eq  27,  are  only  functions  of  x,  and  do  not  affect  the  y- 
direction.  The  reason  for  this,  is  that  the  freestream  flow 
is  parallel  to  the  x-axis.  Because  of  this,  the  y-direction 
contribution  to  the  elemental  upwind  equations  can  be  neglec¬ 
ted.  The  only  term  in  the  elemental  upwind  equations,  given 
by  eqs  28  and  29,  that  is  related  to  the  y-direction  is  the 
BU^j  term.  Whenever  this  term  was  included  in  the  elemental 
stiffness  matrix,  the  solution  diverged,  or  at  best  oscillated 
about  some  solution. 


Test  Coefficient  Influence.  The  test  coefficients  a, 

8  appear  in  eq  25  as 

Wi  =  Nv  +  ttUi  for  hyperbolic  elements 

and 


wv=  Ni+/3Uv  for  transition  elements 

Their  purpose  was  to  vary  the  strength  of  the  upwind  function 
in  order  to  achieve  convergence.  The  full  influence  of  a 
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and  8  cannot  be  expressed  without  considering  the  upwind 
function  itself,  but  some  generalizations  can  be  made.  First, 
it  was  found  that  whenever  either  of  these  values  was  greater 
than  one  or  less  than  zero,  the  solutions  for  any  upwind 
function  diverged  rapidly.  Also,  whenever  the  value  of  8 
for  transition  elements  was  greater  than  ct  for  hyperbolic 
elements,  the  solution  again  diverged  rapidly.  The  solutions 
that  used  a  value  of  8,  one  half  the  value  of  ct  seemed  to  have 
better  solutions,  but  this  depended  on  the  mach  number  and  the 
upwind  function  used.  These  effects  will  be  incorporated  in 
the  subsection  on  different  upwind  functions. 

Mach  Humber  Effects.  The  value  of  the  freestream  mach 
number  Ma,  had  a  large  effect  on  the  convergence  of  the  solu¬ 
tion.  When  Moo>0.87  ,  convergence  never  occurred  for  all 
upwind  functions  and  test  coefficients  tried.  For  Mro<0.87, 
convergent  solutions  were  obtained,  but  depended  on  the 
upwind  function  parameters  and  test  coefficients.  The  de¬ 
tails  are  presented  in  the  next  subsection. 

Upwind  Function  Influence.  The  equations  for  the  up¬ 
wind  function,  given  by  eq  27,  are 

^(1-H,f+ k.  =  o  -i<f<o 

k2  =  o  o<T<i 

The  vertex  of  the  piecewise  parabolas  (H,K)  can  vary. 
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Examples  of  different  upwind  function  parameters  are  shown  in 
Figure  22.  For  explanation  purposes,  the  area  above  the  £- 
axis  will  be  called  the  positive  area,  and  the  area  under  the 
£-axis  will  be  referred  to  as  the  negative  area. 

In  all  three  cases  of  Figure  22,  upwinding  is  present, 
whether  it  be  more  influence  from  the  upstream  part  of  an 
element,  or  less  influence  from  the  downstream  part.  It 
was  found,  that  whenever  the  amount  of  positive  area  was 
greater  than  the  amount  of  negative  area,  Figure  22a,  a 
divergent  solution  occurred.  This  happened  no  matter  what 
value  of  a  and  3  were  used.  Pressure  plots  of  divergent 
solutions  are  not  very  informative  except  to  see  where  the 
divergence  took  place;  so  they  will  not  be  shown.  In  most 
cases,  divergence  occurred  only  where  hyperbolic  elements 
were  present. 

When  the  positive  and  negative  areas  were  equal,  Figure 
22b,  the  solution  was  found  to  be  oscillatory  for  certain 
values  of  a  and  B.  For  a  less  than  one,  all  solutions  di¬ 
verged  no  matter  what  value  of  3  was  used.  When  a  was  equal 
to  one,  and  0.1  <  8  <  0.9,  oscillatory  solutions  for  0.25  5 

<  0.5  and  divergent  solutions  for  other  values  where 
found.  The  oscillatory  solution  in  Figure  23  is  representa¬ 
tive  of  what  happened  for  equal  area  upwind  functions  with 
a=l,  0.1  <  6  <  0.9  and  0.25  <  <  0.5. 

When  the  negative  area  was  greater  than  the  positive  area 
Figure  22c,  convergent  solutions  occurred  for  certain  values 
of  Hj  ,u,  3  and  MOT.  Varying  resulted  in 


a)  More  Positive  Influence  on  Upstream  Nodes 


b)  Positive  Influence  on  Upstream  Nodes  with 

Equal  Negative  Influence  on  Downstream  Nodes 


c)  More  Negative  Influence  on  Downwind  Nodes 


Figure  22. 


Upwind  Functions  for  Three  Different 
Parameters  of  ( li , K ) 


Well  ford  FEM  (Ref  3) 


Figure  23. 


Oscillatory  Solution  for  Equal  Area 
Upwind  Function  (Grid  A):  Mcc-0.86 


one  combination  that  converged,  and  is  shown  in  Figure  22c. 
Using  this  combination  converged  solutions  occurred  at  Meo= 
0.86,  o=l  .0,  8=0.5  and  M^O.86,  ot=0.5,  8=0.5.  All  other  com¬ 
binations  diverged  or  oscillated.  These  two  converged  solu¬ 
tions  are  shown  in  Figures  24  and  25. 

Alternative  Integration  Method 

The  alternative  integration  method  discussed  in  Chapter 
III  was  tried  for  a  couple  of  cases.  These  included  varying 
the  mach  number  and  the  y-direction  influence.  The  y-direc- 
tion  influence  was  varied  by  multiplying  the  BU^j  term  of 
eq  29  by  a  coefficient  ranging  from  zero  to  one.  The  mach 
number  was  ranged  from  0.84  to  0.95. 

In  all  cases,  extremely  divergent  solutions  resulted. 

The  supersonic  region  expanded  to  the  far  field  or  physically 
unrealistic  flow  developed  (i.e.  supersonic  flow  developed 
upstream  of  the  airfoil).  Even  when  more  elements  were 
added  to  the  grid  over  the  airfoil,  divergent  solutions 
resulted.  Figures  26  and  27  are  representative  of  the 
effects  of  the  grid  size  on  the  solution.  Similar  solutions 
resulted  when  the  mach  number  and  y-influence  were  varied. 
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FEM  (Grid  B) 


Figure  24.  Converged  Solution  for  M  =0.86,  a=1.0, 
S=0.5 
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Figure  25.  Converged  Solution  for  Mm=0 . 86 ,  a=0.5, 


26.  Divergent  Solution  Using  Alternative 
Integration  Method  for  Grid  A 
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Figure  27. 


Divergent  Solution  Using  Alternative 
Integration  Method  for  Grid  B 


V. 


Conclusions  and  Recommendations 


The  upwind  and  alternative  integration  methods  used  in 
this  study  did  not  produce  acceptable  solutions  to  the 
small-disturbance  potential  equation  for  transonic  flow. 
When  the  supersonic  region  covered  more  than  two  elements, 
the  solution  diverged,  except  for  the  cases  of  upwind 
parameters  (H,K)  shown  in  Figures  24  and  25. 


Upwind  Function 

The  upwind  method 
recommended  for  futur 
variables  that  must  b 
The  solution  depends 
the  upwind  parameters 
coefficients  a  and  8. 
in  the  finite  element 
sible  to  use  differen 
linear'  or  trigonometr 
parabolic,  and  furthe 


,  described  in  Chapter  III,  is  not 
e  use,  because  there  are  too  many 
e  optimized  to  get  converged  solutions 
on  the  freestream  mach  number  M  , 
(Hj,,K^),  (F^,*^),  and  the  test 
All  of  these  must  be  incorporated 
assembly  process.  It  might  be  pos- 
t  upwind  functions,  such  as  piecewise 
ic  functions  instead  of  piecewise 
r  study  is  warranted  in  this  area. 


Alternative  Integration  Method 

This  method  gave  divergent  results  for  all  cases,  and 
as  it  stands,  is  not  recommended  for  further  study.  Other 
approaches  might  be  taken,  for  example,  assembling  the 
global  stiffness  matrix  an  equation  at  a  time  instead  of  an 
element  at  a  time,  the  effect  of  nodes  that  have  no  in¬ 
fluence  on  the  solution  could  be  zeroed  out.  Another 
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approach,  could  be  to  assume  that  the  border  between  the 
super  and  subsonic  regions  be  considered  a  boundary.  Then 
solve  the  elliptical  part,  and  using  the  solution  on  the 
border,  apply  it  as  a  boundary  condition  when  solving  the 
hyperbolic  part.  This  process  would  have  to  be  done  itera 
tively  until  the  solution  agreed  on  the  boundary  between 
the  regions. 


r - 1 
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Appendix  A 

Finite  Element  Equations  for  a 
Bil  inear  Rectangular  Element 

Elemental  Equations 

The  finite  element  equation  for  the  governing  differen¬ 
tial  equation  in  elemental  form  as  derived  in  Chapter  III 
is 

+  M*Dij 

+M^E,<r>]<tr  -  F-,  <A_1> 

where  A  —  ,  ,  Cij  >  Dij  »  Eij  and  Fi  are  defined  by  eqs  24. 

Shape  Functions 

The  shape  functions  (Ref  1)  for  the  bilinear  rectangular 
element  shown  in  Figure  7  are  given  by 

Ni=i(i+Tit)0+M)  i  =  i,a,v-  <A'2> 

where  ( ? ^ , n £ )  are  the  coordinates  of  the  corner  nodes  in  a 
local  coordinate  system  (£,ri).  Since  the  elemental  equations 
are  expressed  in  global  coordinates,  they  must  be  transformed 
to  local  coordinates  by  the  transformation  equations 

f  =  ^  and  (A-3) 
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Substituting  eqs  A-3  into  the  elemental  equations  given  by 
eq2  24  result  in 

i  i 

A<i=ir//Ni.T  dTdI? 

-i 


B‘i=bjjNw?Ni.7dfd? 

H  -1 

i  I 

cij(4>v  jp  jj Ni,t  Ni>t  NK>  ^  df  dr? 

"l  “1  (A-4) 

I 

t>li=/Nj,TNlf„(t)|d1 

-I 

» 

'*  n=‘« 

i 

p.=  a  f  N,  J„tt)  |  a  T 
A  *i*-\ 


where  j ,  ,  and  apply  only  to  elements  on  the  airfoil. 

In  eqs  A-4,  the  A^j  ,  B—  and  terms  depend  only  on 

the  shape  of  the  elements  and  not  on  their  position.  They 
can  be  integrated  and  evaluated  for  i  =  1,2, 3,4  and  j  =  1,2 
3,4  to  produce  the  following  symmetric  matrices: 

r  -» 
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where 


1  % 

-i 

-  a 

1 

-1 

~  3b 

1 

h 

1 

C11  "Cll 

-C14 

C14 

C^(f)= 

C11 

C14 

-C14 

-r 

u33 

°33 

- 

C33 

c"  =  a?  [  -  C  + 

i<« 

V 

»  r  m  in  i 

in  ii 

C»4~24a*  IT,  “  %  +  1 

tv-4 

,1 

CS3"  "“4^ 

-C] 

where  the  4> '  s  are  the  solutions  at  node  i  from  the  previous 
iteration . 

The  remaining  equations  in  A-4  are  dependent  on  the 
position  of  the  element.  The  quantities  >  Ejj  and  f^ 
exist  only  for  elements  that  border  the  airfoil.  These 
quantities  are  evaluated  at  n  =  “1,  because  a  thin  airfoil 
is  being  used.  The  equation  of  the  airfoil  is 

Y=fM*r<i-4xz)  ;  |xu£ 
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Appendix  B 


Upwind  Functions  and  Elemental  Equations 


In  the  upwind  function  method  described  in  the  text, 
the  weight  function  is  replaced  by 

Ni=  Nt  +  ocuv  (B-l) 

in  elemental  equation  23,  forming  eq  28.  The  terms  in 
the  elemental  upwind  stiffness  matrix  and  force  vector 
(eq  29)  with  the  limits  of  integration  are 

«  i 

auWK,^,,^ 

-( -4. 

CU^ffC 

-V  -t 

(B-2) 

i 

-DUij  =  J  Njf<  U\,  f)X  |  dx 

-y  tjso 

t 

EUiiCf)=[^>K>KNjo<U4^\  <>x 

*5*0 

i 

FUi=  I  T  !,]<)* 

*»  «j=o 
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"Trans forming  ecjs  B— 2  to  local  coordinates  and  substituting 


f  -\<T<o 

Ui.=  ] 

I  ^T-hi)+k2](i+ 9,9)  o<r<t 

and 

N^O-WO-W) 

results  in 
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where 


0 

0 


er 


0 

0 


0  0 
0  0 


0 

0 


a2=t 


w 

3 


A3*4 


x 

4 


2W 

3 


JZ 

2 


Q .  —  — _  4.  A  _L  i_  _Jk  4.  _  \tL  4.  S  ,  «r 

°1  ,2  +  4  +  4  ^  16  +  8  >  2+T 

R^- w  +  ,iL_  -5.  +  T 

B2“  6  12  4  4  16  8  2  ' 


n,-X+ Z.  +  X- JL  aX  _  U-_  A  +  t 
b3  12  4  +  4  >6  6  8  2 
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u4-  16  12  4  4  8  6  2 

ci=  (<k-4)(  -f  +  4^) + (i-4^  W  +  V#) 
c2=(<t>l-4K^r +  T#)+ci-4s)(-Hf +  ^) 

°3=  - •^)+(<k-4)('1^  -  ■§■) 


_  Sa  ya  wa  ua  xyt  _7.yt_'JXc  .  v -r 
™  I  ”  ~£*  ~  16  6  8  *2.  4  4  +At' 


D2  = 


ya 

ie 


+ 


wa  _  u a  xxc  ‘Z-Xt.  vx<. .  S3i  .  -r- 

6  8  "  ie"4~4  2  At' 


where 


X  = 


p. 


70 


Appendix  C 

Elemental  Equations  for  the  Alternative 
Integration  Method 


The  elemental  stiffness  matrix  using  the  alternative 
integration  method  is  evaluated  for  hyperbolic  and  trans¬ 
ition  elements  by  integrating  eqs  24  over  the  area  that 
influences  each  node.  This  area  is  the  area  inside  the 
forward  mach  cone  at  position  j.  For  example,  look  at  the 
elemental  stiffness  matrix  terms 


Kij 


11 

K12 

K13 

K14 

21 

CM 

CM 

K23 

K24 

31 

5* 

U> 

N> 

K33 

K34 

41 

K42 

K43 

K44 

for  the  bilinear  rectangular  element.  K^j  is  the  force  at 

i  due  to  a  unit  potential  at  j.  From  this  statement,  the 

influence  of  each  term  in  K..  can  be  evaluated.  For  the 

ij 

element  in  Figure  C~l,  with  mach  cones  shown,  Figures  C-2 
and  C~3  show  which  K^j  terms  are  influenced  by  the  mach 
cones.  Figure  C-2  is  for  hyperbolic  elements  and  Figure 
C-3  is  for  transition  elements. 


Hyperbolic  Element  Stiffness  Matrix 

The  form  of  the  stiffness  matrix  changes  depending  on 
whether  the  slope  of  the  mach  cone  lines  include  one  or  two 
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Influence  of  Alternative  Integration 
Method  on  the  Elemental  Stiffness  Matrix 
!‘>r  Hyp'-rbi'l  i  c  Elements 
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nodes  (see  Figure  C-2).  For  the  case  where  the  slope  c  of 
the  characteristic  line  is  less  than  the  aspect  ratio  a/b 
of  the  element,  the  terms  of  eq  23  when  integrated  over  the 
area  of  corresponding  mach  cones  result  in 
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[(<t>-4)(c3+c‘>)  + 


A„  =  y  -B,  -B?  -  i-B?  -  y  M?  -  j-  B,  M? 
A44=  y  +B2+  Ba  +  yBa  -  y  M|  +  y  BaMa 
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B„Ba,M,,Ma  FROM  EQS  33 
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When  c  is  greater  than  or  equal  to  a/b  the  resulting  stiff¬ 
ness  matrix  terms  of  eq  23  take  the  form 


0 
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where 


r  -  2B,  +  2B*  2B?  _2_  2B,  2 

M,  M*  3M,  3M,2,  "  3M,5  3 

u _  2B(  +  2B,  ,  2B,3  2 

M,  3M,3  3M,3  3 

J=  2B2  2B2  2Bl  |  2 

Mz  3MJ  3M  |  3 

r  2B>  2  2Bg  2Ba  2B*  2 

"  M2  3M2  M*  3Mj.  3 Mi  3 

In  both  of  the  above  cases,  D^j  and  E^-  from  eq  23  are  the 
same  as  derived  in  Appendix  A  with  E^g  =  E^g  =  D^g  =  D^g  =  0.0 

Transition  Element  Stiffness  Matrix 

.  2 
In  the  transition  element,  the  coefficient  i/c  of  eq  30 

changes  from  negative  to  positive.  A  linear  interpolation 

scheme  must  be  employed  in  order  to  find  the  point  n*  where 

the  coefficient  is  aero.  Using 

2C, 

c,-c* 

to  find  this  position,  where  and  C2  specify  the  value  of 
the  coefficient  of  eq  31  at  the  top  3nd  bottom  of  the  element 
respectively.  For  terms  of  the  stiffness  matrix  that  deal 
with  the  rnach  cone  from  node  4,  it  is  possible  that  the 
intercept  at  q=-l  be  less  than  n*.  For  this  case,  a  new 
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value  of  intercept  must  be  used  in  the  limits  of  integration 
This  new  term  n(-l)  is  calculated  from  the  characteristic 
equations  for  n=-l.  Therefore,  for  stiffness  matrix  terms 
that  deal  with  the  mach  cone  at  node  4  will  have  an  upper 
limit  of  integration  given  by 

17**  =  MIN(«7*  ,q(-0) 

Now,  integrating  eqs  23  over  the  corresponding  areas  in 
Figure  C-3  results  in 
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The  D^j  and  E^j  parts  of  the  stiffness  matrix  are  the  same 


as  shown  in  Appendix  B  with  =  E^  =  =  E 
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